package scu.maqiang.elasticity;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class Elasticity3DT42 {
	public Elasticity3DT42(Mesh3T4 mesh) {
		this.mesh = mesh;
		fs = new FES3T42(mesh);
		K = new SRMatrix(fs.GetNdof());
		x = new double[fs.GetNdof()];
		RHS = new double[fs.GetNdof()];
		disp = new double[3][mesh.getNv()];
		strain = new double[6][mesh.getNv()];
		stress = new double[6][mesh.getNv()];
	}
	
	public static void main(String[] args) {
		long startTime = System.nanoTime();
		Mesh3T4 mesh3 = new Mesh3T4().readMesh("mesh3D\\beam2.mesh", MeshFormat.mesh);
		
		Elasticity3DT42 e = new Elasticity3DT42(mesh3);
		e.elastoStaticProblem(BVPType.COMMON, new double[] { 6.8e10, 0.34 }, new double[]{0.0, 0.0, -1000.0},
				             new int[][] {new int[] {3}}, new Direct[] {Direct.All}, new double[][] {new double[] {0.0, 0.0, 0.0}}, 
				             new int[] {1}, new double[] {0.0, 0.0, 0.0}, "MyElasticity3D.dat");
		//mesh3.toTecplot("MeshBeam.dat");
		// mesh3.printVertices();
//		FES3T42 fs = new FES3T42(mesh3);
//		SRMatrix A = new SRMatrix(fs.GetNdof());
//		fs.assembleElasticStiff(new double[] { 6.8e10, 0.34 }, BVPType.COMMON, A);
//		//double[] XX = new double[fs.GetNdof()];
//		//Arrays.fill(XX, 1.0);
//		//System.out.println(MVO.sum(A.mul(XX)));
//		//System.out.println(A);
//		System.out.println("Number Of DOF " + fs.GetNdof());
//		System.out.println("Size Of A = " + A.getNNZ());
//		double[] RHS = new double[fs.GetNdof()];
//		fs.assembleSurfaceForce(new double[] { 0, 0, -1000.0 }, BVPType.COMMON, RHS, 1);
//		System.out.println(MVO.sum(RHS));
//		//System.out.println(MVO.min(RHS) + "\t" + MVO.max(RHS));
//		fs.applyBC_MBN(A, Direct.All, 3);
//		fs.applyBC_MBN(RHS, Direct.All, 0.0, 3);
////		System.out.println("Size Of A = " + K.getNNZ());
//		
////		fs.applyBCZero(A, RHS, Direct.All, BCType.MBN, 3);
//		System.out.println("Size Of A = " + A.getNNZ());
//		System.out.println(MVO.sum(RHS));
//		double[] x = new double[fs.GetNdof()];
//
//		// IterSSolver solver = new IterSSolver(1e-10, 700);
//		// solver.PCGSSOR(A, RHS, x, 1.7, 1);
//		// solver.CG(A, RHS, x);
//
//		NewIterSSolver newSolver = new NewIterSSolver(A);
//		newSolver.PCGSSOR(RHS, x, 1.7, 1);
//		System.out.println("Solution :");
//		System.out.println("min: " + MVO.min(x) + "\tmax: " + MVO.max(x));
//
//		double[][] uvw = new double[3][mesh3.getNv()];
//		fs.extractComponent(x, uvw[0], uvw[1], uvw[2]);
//
//		mesh3.toTecplot("Elasticity3D.dat", uvw);
//		System.out.println("u min: " + MVO.min(uvw[0]) + "\tmax: " + MVO.max(uvw[0]));
//		System.out.println("v min: " + MVO.min(uvw[1]) + "\tmax: " + MVO.max(uvw[1]));
//		System.out.println("w min: " + MVO.min(uvw[2]) + "\tmax: " + MVO.max(uvw[2]));
		long endTime = System.nanoTime();
		System.out.println("Elapse Time: " + (endTime - startTime) / 1e9 + " s");
	}

	public void elastoStaticProblem(BVPType tp, ScalarFunc[] elasticCoefFunc, double[] paramElasticCoef,
			ScalarFunc[] volumeForce, double[] paramVolumeForce, int[][] fixedDispLabel, Direct[] direct,
			ScalarFunc[][] fixedDisp, double[][] paramFixedDisp, int[] tractionLabel, ScalarFunc[] traction,
			double[] paramTraction, String reportFile) {
		fs.assembleStiff(elasticCoefFunc, paramElasticCoef, tp, K);
		if (volumeForce != null) {
			fs.assembleSource(volumeForce, paramVolumeForce, tp, RHS);
		}
		// System.out.println(MVO.toString(RHS));

		if (tractionLabel != null) {
			fs.assembleSurfaceForce(traction, paramTraction, tp, RHS, tractionLabel);
		}

		int n = direct.length;
		ParamCheck.checkEqual(fixedDispLabel.length, direct.length);
		for (int i = 0; i < n; i++) {
			fs.applyBC_MBN(K, direct[i], fixedDispLabel[i]);
		}

		if (fixedDisp != null) {
			for (int i = 0; i < n; i++) {
				fs.applyBC_MBN(RHS, direct[i], fixedDisp[i], paramFixedDisp[i], fixedDispLabel[i]);
			}
		} else {
			for (int i = 0; i < n; i++) {
				fs.applyBC_MBN(RHS, direct[i], 0.0, fixedDispLabel[i]);
			}
		}

		// System.out.println(MVO.toString(RHS));
		solveAndComputeStrainStressAndExport(reportFile);

	}

	public void elastoStaticProblem(BVPType tp, double[] elasticCoef, double[] volumeForce, int[][] fixedDispLabel,
			Direct[] direct, double[][] fixedDisp, int[] tractionLabel, double[] traction, String reportFile) {
		fs.assembleStiff(elasticCoef, tp, K);
		double[] XX = new double[fs.GetNdof()];
		Arrays.fill(XX, 1.0);
		System.out.println(MVO.sum(K.mul(XX)));
		//fs.assembleElasticStiff(new double[] { 6.8e10, 0.34 }, BVPType.COMMON, K);
		if (volumeForce != null) {
			fs.assembleSource(volumeForce, tp, RHS);
		}
		//System.out.println(MVO.toString(RHS));

		if (tractionLabel != null) {
			fs.assembleSurfaceForce(traction, tp, RHS, tractionLabel);
		}
		
		//System.out.println(MVO.sum(RHS));
		//System.out.println(MVO.min(RHS) + "\t" + MVO.max(RHS));
		
		int n = direct.length;
		ParamCheck.checkEqual(fixedDispLabel.length, direct.length);
		for (int i = 0; i < n; i++) {
			fs.applyBC_MBN(K, direct[i], fixedDispLabel[i]);
		}
		
		//System.out.println("Size Of A = " + K.getNNZ());

		if (fixedDisp != null) {
			for (int i = 0; i < n; i++) {
				fs.applyBC_MBN_Disp(RHS, direct[i], fixedDisp[i], fixedDispLabel[i]);
			}
		} else {
			for (int i = 0; i < n; i++) {
				fs.applyBC_MBN(RHS, direct[i], 0.0, fixedDispLabel[i]);
			}
		}
		System.out.println(MVO.sum(RHS));
//		fs.applyBCZero(K, RHS, Direct.All, BCType.MBN, 3);

		// System.out.println(MVO.toString(RHS));
		NewIterSSolver solver = new NewIterSSolver(K);
		solver.PCGSSOR(RHS, x, 1.7, 1);
		fs.extractComponent(x, disp);
		mesh.toTecplot(reportFile, disp);
		double[][] uvw = new double[3][mesh.getNv()];
		fs.extractComponent(x, uvw);

		//mesh3.toTecplot("Elasticity3D.dat", uvw);
		System.out.println("u min: " + MVO.min(uvw[0]) + "\tmax: " + MVO.max(uvw[0]));
		System.out.println("v min: " + MVO.min(uvw[1]) + "\tmax: " + MVO.max(uvw[1]));
		System.out.println("w min: " + MVO.min(uvw[2]) + "\tmax: " + MVO.max(uvw[2]));
	}

	public void solveAndComputeStrainStressAndExport(String reporFile) {
		NewIterSSolver solver = new NewIterSSolver(K);
		solver.PCGSSSOR(RHS, x, 1.7, 1);
		fs.extractComponent(x, disp);
		mesh.toTecplot(reporFile, disp);
		double[][] uvw = new double[3][mesh.getNv()];
		fs.extractComponent(x, uvw);

		//mesh3.toTecplot("Elasticity3D.dat", uvw);
		System.out.println("u min: " + MVO.min(uvw[0]) + "\tmax: " + MVO.max(uvw[0]));
		System.out.println("v min: " + MVO.min(uvw[1]) + "\tmax: " + MVO.max(uvw[1]));
		System.out.println("w min: " + MVO.min(uvw[2]) + "\tmax: " + MVO.max(uvw[2]));
		// fs.computeNodalStrain(x, strain);
		// fs.computeNodalStress(x, stress);
		// fs.computeEqStress(x, eqStress);
		// fs.computeNodalStrainStress(x, strain,stress, eqStress);
	}

	public double[][] exportDisplacement() {
		double[][] myDisp = new double[3][mesh.getNv()];
		myDisp[0] = disp[0].clone();
		myDisp[1] = disp[1].clone();
		myDisp[2] = disp[2].clone();
		return myDisp;
	}

	Mesh3T4 mesh;
	FES3T42 fs;
	SRMatrix K;
	double[] RHS;
	double[] x;
	double[][] disp;
	double[][] strain;
	double[][] stress;
	double[][] eqStress;

}
